

t1a_lm = lm_robust(attitudinal_outcome ~ I(treatment == "T1"), 
                   data = survey_df %>% filter(treatment %in% c("T1", "C"))) %>% tidy()
t1b_lm = lm_robust(behavioral_outcome ~ I(treatment == "T1"), 
                   data = survey_df %>% filter(treatment %in% c("T1", "C"))) %>% tidy()

t2a_lm = lm_robust(attitudinal_outcome ~ I(treatment == "T2"),
                   data = survey_df %>% filter(treatment %in% c("T2", "C"))) %>% tidy()
t2b_lm = lm_robust(behavioral_outcome ~ I(treatment == "T2"), 
                   data = survey_df %>% filter(treatment %in% c("T2", "C"))) %>% tidy()


t3a_lm = lm_robust(attitudinal_outcome ~ I(treatment == "T3"), 
                   data = survey_df %>% filter(treatment %in% c("T3", "C"))) %>% tidy()
t3b_lm = lm_robust(behavioral_outcome ~ I(treatment == "T3"), 
                   data = survey_df %>% filter(treatment %in% c("T3", "C"))) %>% tidy()

t4a_lm = lm_robust(attitudinal_outcome ~ I(treatment == "T4"), 
                   data = survey_df %>% filter(treatment %in% c("T4", "C"))) %>% tidy()
t4b_lm = lm_robust(behavioral_outcome ~ I(treatment == "T4"), 
                   data = survey_df %>% filter(treatment %in% c("T4", "C"))) %>% tidy()


individual_mods = list(t1a_lm, t1b_lm, t2a_lm, t2b_lm,t3a_lm, t3b_lm, t4a_lm, t4b_lm)

individual_hypotheses = 
  bind_rows(individual_mods) %>%
  mutate(val = str_detect(term, "treat")) %>%
  filter(val == T) %>%
  mutate(outcome_lab = case_when(str_detect(outcome, "att") ~ "attitudinal compliance",
                                 str_detect(outcome, "beh") ~ "behavioral compliance")) %>%
  mutate(treatment_lab = case_when(str_detect(term, "1") ~ "T1: Religious Institution",
                                   str_detect(term, "2") ~ "T2: Secular Institution",
                                   str_detect(term, "3") ~ "T3: Secular Personal",
                                   str_detect(term, "4") ~ "T4: Religious Personal")) %>%
  mutate(treatment_lab = factor(treatment_lab, levels = c("T4: Religious Personal", "T3: Secular Personal", 
                                                          "T2: Secular Institution", "T1: Religious Institution")))


ggplot(data = individual_hypotheses) + 
  geom_hline(yintercept=0, lty=2) + 
  geom_pointrange(aes(x=treatment_lab, y=estimate, ymin = (estimate - 1.96*std.error), ymax=(estimate + 1.96*std.error))) + 
  geom_pointrange(aes(x=treatment_lab, y=estimate, ymin = (estimate - 1.65*std.error), ymax=(estimate + 1.65*std.error)), fatten=2, size=1.1) + 
  labs(x='', y = '') + 
  ylab(expression("Average Treatment Effect (" ~ beta[1] ~ ")")) +
  facet_wrap(outcome_lab~.) + 
  coord_flip() + 
  theme_bw() +
  theme(legend.title=element_blank(), text=element_text(size=20), 
        axis.title.y=element_blank()) +
  theme(text=element_text(size=12, family="Times")) + 
  theme(panel.grid.minor = element_blank(), panel.grid.major.x = element_blank(), panel.grid.major.y = element_line(size = .2)) + 
  theme(strip.background =element_rect(fill="white"))


ggsave("./outputs/sa_figure_2.pdf", width = 6, height = 3)

rm(list=setdiff(ls(), "survey_df"))